********************************
* This file constructs Tables 1-4
* Instead of weighting, I control for X
********************************
clear all
set mem 4000m
set matsize 2000
set maxvar 8000

use build.dta


********************************
* NEW VARIABLES
********************************

* Per capita wage in manufacturing, in real dollars
* Manufacturing payroll of production workers over average number of production workers 
g wage0    = (pcmwage00)/(cpi0/100)
g wage20   = (pcmwage20)/(cpi20/100)
g wage30   = (pcmwage30)/(cpi30/100)
g wage40   = (pcmwage39)/(cpi40/100)
g wage50   = (pcmwage47)/(cpi50/100)
g wage60   = (pcmwage58)/(cpi60/100)
g wage70   = (pcmwage67)/(cpi70/100)
g wage80   = (pcmwage82)/(cpi80/100)
g wage90   = (pcmwage87)/(cpi90/100)
g wage2000 = (pcmwage97)/(cpi2000/100)


* Per capita wage in trade, in real dollars
* total payroll in wholesale establishments + retail establishments / (total workers in retail+ in wholsale)
g twage30  = pctwage30/(cpi30/100) 
g twage40  = pctwage40/(cpi40/100)
g twage50  = pctwage54/(cpi50/100)
g twage60  = pctwage63/(cpi60/100)
g twage70  = pctwage72/(cpi70/100)
g twage80  = pctwage82/(cpi80/100)
g twage90  = pctwage87/(cpi90/100)
g twage2000  = pctwage97/(cpi2000/100)

* Agricultural values

g lnfaval0 = ln(faval900/(cpi0/100))
g lnfaval10 = ln(faval910/(cpi10/100))
g lnfaval20 = ln(faval920/(cpi20/100))
g lnfaval30 = ln(faval930/(cpi30/100))
g lnfaval40 = ln(faval940/(cpi40/100))
g lnfaval50 = ln(faval950/(cpi50/100))
g lnfaval60 = ln(faval959/(cpi60/100))
g lnfaval70 = ln(faval1969/(cpi70/100))
g lnfaval80 = ln(faval1982/(cpi80/100))
g lnfaval90 = ln(faval1992/(cpi90/100))
g lnfaval2000=ln(faval2002/(cpi2000/100))


* Median family income
gen lnmedfaminc50  = ln(medfaminc50/(cpi50/100)) 
gen lnmedfaminc60  = ln(medfaminc60/(cpi60/100)) 
gen lnmedfaminc70  = ln(medfaminc70/(cpi70/100)) 
gen lnmedfaminc80  = ln(medfaminc80/(cpi80/100)) 
gen lnmedfaminc90  = ln(medfaminc80/(cpi90/100)) 
gen lnmedfaminc2000= ln(medfaminc2000/(cpi2000/100))

* Farm production
gen lnvfprod30   = log(vfprod30/(cpi30/100))
gen lnvfprod40   = log(vfprod40/(cpi40/100))
gen lnvfprod50   = log(vfprod50/(cpi50/100))
gen lnvfprod60   = log(vfprod60/(cpi60/100))
gen lnvfprod70   = log(vfprod70/(cpi70/100))
gen lnvfprod80   = log(vfprod80/(cpi80/100))
gen lnvfprod90   = log(vfprod90/(cpi90/100))
gen lnvfprod2000 = log(vfprod2000/(cpi2000/100))

*Foreign born
gen fb0=fbwmtot00 + fbwftot00 
gen fb10=fbwtot10
gen fb20=fbwmtot20 + fbwftot20 
gen fb30=fbwmtot + fbwftot

gen fbshr20=fb20/(wmtot20 + wftot20)
gen fbshr30=fb30/(wmtot + wftot)

*Housing Values/Rents

ren medrnt30_NHGIS medrnt30
ren medhsval30_NHGIS medhsval30
ren var88_county72 medhsval70
ren var89_county72 medrnt70

foreach var in medhsval medrnt{
	foreach yr in 30 40 50 60 70 80 90 2000{
		cap replace `var'`yr'=`var'`yr'/(cpi`yr'/100)
	}
}

*various
drop other60

*drop counties experiencing big changes in area
gen d=(b1_lnd01_county00 - area)/(b1_lnd01_county00 + area)/2
drop if abs(d)>.03
replace area=(b1_lnd01_county00 + area)/2



******************************
* standardize variable names *
******************************


ren emp00 emp0
ren manuf_jobs_00 manuf_jobs_0

foreach yr in 0 10 20 30 40 50 60 70 80 90 2000{
	cap drop manuf`yr'
	cap drop agr`yr'
	ren manuf_jobs_`yr' manuf`yr'
	cap ren ag_jobs`yr' agr`yr'
}




foreach yr in 0 10 20 30 60 80 90 2000{
	gen other`yr'=emp`yr'-agr`yr'-manuf`yr'
}


*********************
*    Make Share     *
*********************

foreach var in manuf agr{
	foreach yr in 0 10 20 30 40 50 60 70 80 90 2000{
		cap gen `var'shr`yr'=`var'`yr'/emp`yr'
	}
}

********************
* prepare outcomes *
********************


foreach var in pop emp house wage twage agr manuf other medhsval medrnt fb{
	cap gen ln`var'0=ln(`var'0)
	cap gen ln`var'10=ln(`var'10)
	cap gen ln`var'20=ln(`var'20)
	cap gen ln`var'30=ln(`var'30)
	cap gen ln`var'40=ln(`var'40)
	cap gen ln`var'50=ln(`var'50)
	cap gen ln`var'60=ln(`var'60)
	cap gen ln`var'70=ln(`var'70)
	cap gen ln`var'80=ln(`var'80)
	cap gen ln`var'90=ln(`var'90)
	cap gen ln`var'2000=ln(`var'2000)
}

drop lntwage20  lnmedhsval20 lnmedrnt20 //stata confuses 20 and 2000


*******************
*Generate Controls*
*******************
g urbshare0=popurb0/pop0
g urbshare10=popurb10/pop10
g urbshare20=popurb20/pop20
g urbshare30=popurb30/pop20
gen popdens0=pop0/b1_lnd01_county00


*fix zeros in covariate quantities

foreach var in agr manuf other{
	foreach yr in 10 20 30{
		cap replace ln`var'`yr'=0 if `var'`yr'<=0
		gen no`var'`yr'dum=`var'`yr'<=0
	}
}

ren mfgcap00 mfgcap0
gen lnmfgcap0=ln(mfgcap0)
replace lnmfgcap0=0 if mfgcap0==0



*fix wages

foreach var in wage twage{
	foreach yr in 20 30{
		cap replace ln`var'`yr'=-1 if `var'`yr'==.
		cap gen no`var'`yr'dum=`var'`yr'==.
	}
}



*transformations
foreach var of varlist elevmax elevrang popdens0 tillit10 tillit1020 tillit1010 retsales radiorep totunemp area{
	gen ln`var'=ln(`var')
}


center tmean* lnelevmax

foreach var of varlist c_lnelevmax white0 white20 white30 c_tmean* lnmanuf0 lnmanuf10 lnmanuf20 lnradiorep lnemp20 lnemp30 lnwage0 lnwage20 lnwage30 lntwage30 lnpop0 lnpop20 lnpop30 lnfaval0 lnfaval20 lnfaval30 tmean*{
	gen `var'sq=`var'^2
	gen `var'cub=`var'^3
}

gen popdifsq=(lnpop30-lnpop20)^2
gen empdifsq=(lnemp30-lnemp20)^2
gen manufdifsq=(lnmanuf20-lnmanuf0)^2
gen urbsharedifsq=(urbshare20-urbshare0)^2
gen whitedifsq=(white20-white0)^2
gen agrdifsq=(lnagr20-lnagr0)^2
gen wagedifsq=(lnwage20-lnwage0)^2
gen favaldifsq=(lnfaval20-lnfaval0)^2

gen agrshr30sq=agrshr30^2
gen agrshr20sq=agrshr20^2

gen lnagr20sq=lnagr20^2
gen lnagr0sq=lnagr0^2
gen fbdifsq=(lnfb20-lnfb0)^2

gen wagedif2sq=(lnwage30-lnwage20)^2
gen tmean_jan_jul=tmean_jan*tmean_jul

gen pctil20=tillit1020/t10tot20
gen pctil30=tillit10/t10tot
replace PRADIO=PRADIO/100
gen urate30=totunemp/(totunemp+emp30)



global X "lnelevmax lnelevrang lnarea lnpop20 lnpop20sq lnpop30 lnpop30sq popdifsq agrshr20 agrshr20sq agrshr30 agrshr30sq manufshr20 manufshr30 nowage20 nowage30 lnwage20 lnwage30 notwage30 lntwage30 lnemp20 lnemp30 urbshare20 urbshare30 lnfaval20 lnfaval30 lnmedhsval30 lnmedrnt30 white20 white20sq white30 white30sq pctil20 pctil30 PRADIO urate30 fbshr20 fbshr30"

run "~/geo/x_ols_JP_v11.ado" //use Juan Pablo's version of Conley (1999)
global tables="~/latex/tables_raw"


*******************************************
*******************************************
* Propensity score 
*******************************************
*******************************************

sum $X

logit tva $X, cluster(state)
predict phat

keep if e(sample)==1

sum phat if tva==1, det
sum phat if tva==0, det
local cut=r(p25)

tab tva

preserve
drop if phat<`cut'&tva==0

**Construct O-B weights**

mata: D=st_data(.,"tva")
mata: one=J(rows(D),1,1)
mata: nD=one-D
order $X
local K:word count $X
mata: X=st_data(.,1..`K')
mata: X=(one, X)
mata: w=D'*X*invsym(quadcross(X,nD,X))*X'/sum(D)
mata: w=w':*nD
gen w=.
mata: st_store(.,"w",w)

replace w=1 if tva==1
sum w if tva==0
sum w if tva==1


/*
**Make Maps**
lab var phat "Probability of Treatment"
lab var tva "TVA"
spmap w using county_coordinates if tva==0&phat>`cut'&phat<., id(id) title("Control Counties in Estimation Sample") legtitle("Oaxaca Weight")
gr export /tmp/map0.eps, replace
spmap w using county_coordinates if tva==1&phat!=., id(id) title("TVA Counties in Estimation Sample")
gr export /tmp/map1.eps, replace
!ps2pdf -dEPSCrop /tmp/map0.eps map0.pdf
!ps2pdf -dEPSCrop /tmp/map1.eps map1.pdf
!rm /tmp/map0.eps /tmp/map1.eps

spmap tva using county_coordinates if phat>`cut'&phat<., id(id) title("Counties in Estimation Sample") legend(off) note("Shaded Areas are in TVA")
gr export /tmp/fullmap.eps, replace
!ps2pdf -dEPSCrop /tmp/fullmap.eps fullmap.pdf
!rm /tmp/fullmap.eps
*/
restore





***************************
** Table 1 Summary Stats 
***************************


gen lnwageobs30=ln(wage30)
gen lnwageobs20=ln(wage20)
foreach var in lnpop lnemp lnhouse lnwageobs manufshr agrshr white urbshare pctil fbshr lnfaval {
	gen D_`var'=`var'30-`var'20
	lab var D_`var' "Change in `var' 1920-1930"
}


qui estpost tabstat lnpop30 lnemp30 lnhouse30 lnwageobs30 manufshr30 agrshr30 white30 urbshare30 pctil30 fbshr30 lnfaval30 lnmedhsval30 lnmedrnt30 PRADIO elevmax elevrang south D_* if tva==1, stat(mean) col(stat)
eststo sum0
qui inspect state if tva==1
estadd scalar clusters = r(N_unique)


qui estpost tabstat lnpop30 lnemp30 lnhouse30 lnwageobs30 manufshr30 agrshr30 white30 urbshare30 pctil30 fbshr30 lnfaval30 lnmedhsval30 lnmedrnt30 PRADIO elevmax elevrang south D_* if tva==0, stat(mean) col(stat)
eststo sum1 
qui inspect state if tva==0
estadd scalar clusters = r(N_unique)

qui estpost tabstat lnpop30 lnemp30 lnhouse30 lnwageobs30 manufshr30 agrshr30 white30 urbshare30 pctil30 fbshr30 lnfaval30 lnmedhsval30 lnmedrnt30 PRADIO elevmax elevrang south D_* if tva==0&south==1, stat(mean) col(stat)
eststo sum2 
qui inspect state if tva==0&south==1
estadd scalar clusters = r(N_unique)


drop if phat<`cut'


qui estpost tabstat lnpop30 lnemp30 lnhouse30 lnwageobs30 manufshr30 agrshr30 white30 urbshare30 pctil30 fbshr30 lnfaval30 lnmedhsval30 lnmedrnt30 PRADIO elevmax elevrang south D_* if tva==1, stat(mean) col(stat)
eststo sum3
qui inspect state if tva==1
estadd scalar clusters = r(N_unique)

qui estpost tabstat lnpop30 lnemp30 lnhouse30 lnwageobs30 manufshr30 agrshr30 white30 urbshare30 pctil30 fbshr30 lnfaval30 lnmedhsval30 lnmedrnt30 PRADIO elevmax elevrang south D_* if tva==0, stat(mean) col(stat)
eststo sum4 
qui inspect state if tva==0
estadd scalar clusters = r(N_unique)

qui estpost tabstat lnpop30 lnemp30 lnhouse30 lnwageobs30 manufshr30 agrshr30 white30 urbshare30 pctil30 fbshr30 lnfaval30 lnmedhsval30 lnmedrnt30 PRADIO elevmax elevrang south D_* if tva==0&south==1, stat(mean) col(stat) 
eststo sum5 
qui inspect state if tva==0&south==1
estadd scalar clusters = r(N_unique)


esttab sum0 sum1 sum2 sum3 sum4 sum5 using $tables/summ.csv, cells(mean(fmt(3))) stat(N clusters, fmt(0)) replace mlabels(lnpop30 lnemp30 lnhouse30 lnwageobs30 manufshr30 agrshr30 white30 urbshare30 pctil30 fbshr30 lnfaval30 lnmedhsval30 lnmedrnt30 PRADIO elevmax elevrang south) label unstack nodepvar

******************************
*****************************
* TABLE 2 -- Selection
*****************************
******************************


local i=1
foreach var in lnpop lnemp lnhouse lnwage manufshr agrshr lnfaval{	
	g  Dy = (`var'40-`var'0)/4
	*Winsorize outliers
	qui centile Dy, c(1 99)
	qui replace Dy=r(c_1) if Dy<r(c_1)
	qui replace Dy=r(c_2) if Dy>r(c_2)&Dy!=.

	qui reg Dy tva, cluster(state) 
	eststo base`i'
	qui x_ols Dy $X, cluster(state) bo(tva)
	eststo eval`i'
	qui x_ols Dy $X, lat(latitude) long(longitud) cut1(200) cut2(200) bo(tva)
	eststo cor_eval`i'
	drop Dy 
	local ++i
}



esttab base* using $tables/base1.csv, replace keep(tva) mlabels(lnpop lnemp lnhouse lnwage manufshr agrshr lnfaval) b(3) se(3) star(* 0.1 ** 0.05 *** 0.01)
esttab eval* using $tables/eval1.csv, replace keep(tva) mlabels(lnpop lnemp lnhouse lnwage manufshr agrshr lnfaval) b(3) se(3) star(* 0.1 ** 0.05 *** 0.01)
esttab cor_* using $tables/cor1.csv, replace keep(tva) mlabels(lnpop lnemp lnhouse lnwage manufshr agrshr lnfaval) b(3) se(3) star(* 0.1 ** 0.05 *** 0.01)





******************************
*****************************
* TABLE 2A -- Selection (South Only)
*****************************
******************************

local i=1
foreach var in lnpop lnemp lnhouse lnwage manufshr agrshr lnfaval{	
	g  Dy = (`var'40-`var'0)/4

	qui centile Dy, c(1 99)
	qui replace Dy=r(c_1) if Dy<r(c_1)
	qui replace Dy=r(c_2) if Dy>r(c_2)&Dy!=.

	qui x_ols Dy tva if south==1, lat(latitude) long(longitud) cut1(200) cut2(200)
	eststo base`i', title(`var')
	qui x_ols Dy $X if south==1, lat(latitude) long(longitud) cut1(200) cut2(200) bo(tva)
	eststo cor_eval`i', title(`var')
	drop Dy 
	local ++i
}



esttab base* using $tables/base1a.csv, replace keep(tva) mlabels(lnpop lnemp lnhouse lnwage manufshr agrshr lnfaval) b(3) se(3) star(* 0.1 ** 0.05 *** 0.01)
esttab cor_* using $tables/cor1a.csv, replace keep(tva) mlabels(lnpop lnemp lnhouse lnwage manufshr agrshr lnfaval) b(3) se(3) star(* 0.1 ** 0.05 *** 0.01)


******************************
*****************************
* TABLE 3 -- Impacts
*****************************
******************************



local i=1
foreach var in lnpop lnwage lnagr lnmanuf lnvfprod lnmedfaminc lnfaval lnmedhsval{
	
	if "`var'"=="lnmedfaminc"{
		g Dy=(`var'2000-`var'50)/5
	}
	else{
		g  Dy = (`var'2000-`var'40)/6
	}

	qui centile Dy, c(1 99)
	qui replace Dy=r(c_1) if Dy<r(c_1)
	qui replace Dy=r(c_2) if Dy>r(c_2)&Dy!=.

	qui reg Dy tva, cluster(state) 
	eststo base`i', title(`var')
	qui x_ols Dy $X, cluster(state) bo(tva)
	eststo eval`i', title(`var')
	qui x_ols Dy $X, lat(latitude) long(longitud) cut1(200) cut2(200) bo(tva)
	eststo cor_eval`i', title(`var')
	drop Dy 
	local ++i
}


esttab base* using $tables/base2.csv, replace keep(tva) mlabels(lnpop lnwage lnagr lnmanuf lnvfprod lnmedfaminc lnfaval lnmedhsval) b(3) se(3) star(* 0.1 ** 0.05 *** 0.01)
esttab eval* using $tables/eval2.csv, replace keep(tva) mlabels(lnpop lnwage lnagr lnmanuf lnvfprod lnmedfaminc lnfaval lnmedhsval) b(3) se(3) star(* 0.1 ** 0.05 *** 0.01)
esttab cor_* using $tables/cor2.csv, replace keep(tva) mlabels(lnpop lnwage lnagr lnmanuf lnvfprod lnmedfaminc lnfaval lnmedhsval) b(3) se(3) star(* 0.1 ** 0.05 *** 0.01)



******************************
***  Table 3A (South Only) ***
******************************

local i=1
foreach var in lnpop lnwage lnagr lnmanuf lnvfprod lnmedfaminc lnfaval lnmedhsval{
	
	if "`var'"=="lnmedfaminc"{
		g Dy=(`var'2000-`var'50)/5
	}
	else{
		g  Dy = (`var'2000-`var'40)/6
	}

	qui centile Dy, c(1 99)
	qui replace Dy=r(c_1) if Dy<r(c_1)
	qui replace Dy=r(c_2) if Dy>r(c_2)&Dy!=.

	qui x_ols Dy tva if south==1, lat(latitude) long(longitud) cut1(200) cut2(200)
	eststo base`i', title(`var')
	qui x_ols Dy $X if south==1, lat(latitude) long(longitud) cut1(200) cut2(200) bo(tva)
	eststo cor_eval`i', title(`var')
	drop Dy 
	local ++i
}


esttab base* using $tables/base2a.csv, replace keep(tva) mlabels(lnpop lnwage lnagr lnmanuf lnvfprod lnmedfaminc lnfaval lnmedhsval) b(3) se(3) star(* 0.1 ** 0.05 *** 0.01)
esttab cor_* using $tables/cor2a.csv, replace keep(tva) mlabels(lnpop lnwage lnagr lnmanuf lnvfprod lnmedfaminc lnfaval lnmedhsval) b(3) se(3) star(* 0.1 ** 0.05 *** 0.01)




*****************************************************
*****************************************************
* TABLE 4
* By period 30-60 60-2000
*****************************************************
*****************************************************

local i=1
foreach var in lnpop lnwage lnagr lnmanuf lnvfprod lnmedfaminc lnfaval lnmedhsval{
	g Dy=(`var'2000-`var'60)/4


	qui centile Dy, c(1 99)
	qui replace Dy=r(c_1) if Dy<r(c_1)
	qui replace Dy=r(c_2) if Dy>r(c_2)&Dy!=.

	qui x_ols Dy $X, cluster(state) bo(tva)
	eststo eval1_`i', title(`var')
	qui x_ols Dy $X if south==1, lat(latitude) long(longitud) cut1(200) cut2(200) bo(tva)

	eststo eval1a_`i', title(`var')

	if "`var'"!="lnmedfaminc"{
		replace Dy=(`var'60-`var'40)/2
		qui centile Dy, c(1 99)
		qui replace Dy=r(c_1) if Dy<r(c_1)
		qui replace Dy=r(c_2) if Dy>r(c_2)&Dy!=.

		qui x_ols Dy $X, cluster(state) bo(tva)
		eststo eval2_`i', title(`var')
		qui x_ols Dy $X if south==1, lat(latitude) long(longitud) cut1(200) cut2(200) bo(tva)
		eststo eval2a_`i', title(`var')
	}

	drop Dy 
	local ++i
}


esttab eval1_* using $tables/eval3first.csv, replace keep(tva) mlabels(lnpop lnwage lnagr lnmanuf lnvfprod lnmedfaminc lnfaval lnmedhsval) b(3) se(3) star(* 0.1 ** 0.05 *** 0.01)
esttab eval1a_* using $tables/eval3firsta.csv, replace keep(tva) mlabels(lnpop lnwage lnagr lnmanuf lnvfprod lnmedfaminc lnfaval lnmedhsval) b(3) se(3) star(* 0.1 ** 0.05 *** 0.01)

esttab eval2_* using $tables/eval3second.csv, replace keep(tva) mlabels(lnpop lnwage lnagr lnmanuf lnvfprod lnmedfaminc lnfaval lnmedhsval) b(3) se(3) star(* 0.1 ** 0.05 *** 0.01)
esttab eval2a_* using $tables/eval3seconda.csv, replace keep(tva) mlabels(lnpop lnwage lnagr lnmanuf lnvfprod lnmedfaminc lnfaval lnmedhsval) b(3) se(3) star(* 0.1 ** 0.05 *** 0.01)


